Periodic orbit theory for realistic cluster potentials: 
The leptodermous expansion 
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The formation of supershells observed in large metal clusters can be qualitatively understood 
from a periodic-orbit-expansion for a spherical cavity. To describe the changes in the supershell 
structure for different materials, one has, however, to go beyond that simple model. We show 
how periodic-orbit-expansions for realistic cluster potentials can be derived by expanding only the 
classical radial action around the limiting case of a spherical potential well. We give analytical 
results for the leptodermous expansion of Woods-Saxon potentials and show that it describes the 
shift of the supershells as the surface of a cluster potential gets softer. As a byproduct of our work, 
we find that the electronic shell and supershell structure is not affected by a lattice contraction, 
which might be present in small clusters. 
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I. INTRODUCTION 



One of the most surprising aspects of the physics of 
metal clusters is the supershell structure observed in 
mass-abundance spectraEfQ This feature can be traced 
back to a beating pattern in the density of states for 
typical cluster potentials!! The conceptual framework 
for understanding how this quantum interference comes 
about is provided by periodic-orbit-theory.&El The ele- 
gance of this approach rests on the fact that the periodic- 
orbit-expansion (POE) is known analytically for the 
spherical cavity. For this model potential it was found 
that the most important contributions to the oscillating 
part of the density of states stem from the two short- 
est planar periodic orbits: the triangular and the square 
orbits. Since these contributions oscillate with similar 
frequencies, their interference gives rise to a beating pat- 
tern, hence supershells. 

Although a spherical potential well is a good first ap- 
proximation to a cluster potential, this model clearly can- 
not account for the changes in the electronic shell and 
supershell structure observed for clusters made of dif- 
ferent materials. It is therefore desirable to understand 
how the periodic-orbit-expansion is modified as one con- 
siders more realistic model potentials. A straightforward 
approach for doing so is to solve the action integrals, 
which lie at the heart of periodic-orbit-theory, numer- 
ically. That way, however, most of the elegance and 
power of the periodic-orbit-expansion is lost. Analytical 
expression, on the other hand, may well reveal the rel- 
evant parameters determining the supershell structure, 
and provide insight into how it changes as the cluster 
potential is varied. They should prove especially help- 
ful in the search for better self-consistent models, which 
properly describe the experimental data. 

While, the simple spherical, homogeneous jellium 
modeH'El works quite well for the alkali clusters, it fai 
to describe the supershell structure observed in G&n 



Attempts to improve the situation iridude, e.g. the in- 
troduction of smpot±r|jellium-profiles,Lj the inclusion of 
pseudopoteiitials,t2rllj or the consideration of surface 
roughness S3 

The present work arises from the desire to understand 
the simple spherical, homogeneous jellium model. An 
analysis of the density dependence of the electronic su- 
pershells in jellium clusters showed that the supershells 
are shifted as-,the potential at the cluster surface be- 
comes softer.Ej It has been demonstrated that this shift 
can be understood in the framework of a periodic-orbit-, 
expansion for typical self-consistent cluster potentials.E£l 
The purpose of the present paper is to give a derivation 
of the leptodermous expansion, which requires the lin- 
earization of only the radial action. We furthermore ana- 
lyze the validity of the approximations involved and show 
comparisons with quantum mechanical calculations. 

To set the stage for the semiclassical treatment of elec- 
tronic supershell^ Sep. |TJ gives a review of the shell cor- 
rection methods.li3~ll3 These methods establish a system- 
atic relation between self-consistent calculations and one- 
electron calculations for suitable model potentials. We 
stress the fact that it is decisive to choose families of po- 
tentials that vary smoothly with cluster size, to describe 
the electronic shell structure properly. 

Sec. [II is devoted to periodic-orbit-expansions. We 



sketch the derivation of the POE for the oscillating part p 
of the density of states using the path integral formalism 
along the lines given by Gutzwiller.0 Special attention is 
paid to the rate of convergence of the sum over classical 
periodic orbits. Since the shell and supershell structure 
observed in the mass spectra of metal clusters is not di- 
rectly linked to p but rather the variations E in total 
energy, we proceed to derive a periodic-orbit-expansion 
for E. We find that the latter expansion converges much 
more rapidly than that for the density of states, hence 
making any artificial smoothing of the spectrum, com- 
monly introduced to lessen the contribution of the longer 
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orbits to pJliS superfluous. To assess the validity of the 
expression for E we check the truncated POE against the 
quantum- mechanical result. This comparison shows that 
in the size range, which seems experimentally accessible, 
E is well described by a truncated POE, taking only tri- 
angular and square orbits into account. This justifies the 
common practice of truncating the POE after the two 
shortest planar orbits. 

As an immediate application of the periodic-orbit- 
expansion for E we show in Sec. |lV] that the electronic 
shell structure is virtually not affected by any lattice con- 
traction. Such an increase in density_£pr small clusters 
was sugg ested by EXAFS analysesBEl and seems also 
plausible from the viewpoint of continuum mechanics, 
which implies that the surface tension should lead to a 
compression of the smaller clusters. We show that the 
major effect of a lattice contraction on E is a change of 
the amplitude for small clusters. The position of the shell 
minima and the supershell structure is, however, not no- 
ticeably changed, even for unrealistically large contrac- 
tions. 

In Sec. [v| we show how to extend the periodic-orbit- 
expansion of E to more realistic potentials. We start 
from the observation that the surface width a of the clus- 
ter potential is an imaortant parameter determining the 
supershell structure. E-3 The basic idea is then to expand 
the action integrals entering the POE around the ana- 
lytically known results for a potential well. It turns out 
that the actions can be very well approximated by lin- 
ear functions in the surface parameter a/ Rq, where Rq is 
the radius of the cluster. Thus a finite surface width 
leaves the frequencies in the POE unchanged and, to 
first order, only introduces phase shifts. Taking also the 
change of the Fermi energy into account, we can under- 
stand the changes in the electronic shells and supershells 
introduced by a soft potential surface. 

The technical details of the leptodermous expansion for 
Woods-Saxon potentials are described in the appendix. 

To simplify the notation we set h 2 /2m to unity, i.e. we 
give lengths in Bohr-radii (do) and energies in Rydberg. 



II. SHELL CORRECTION METHODS 

The total energy E(N) of clusters having N valence 
electrons can be split into a smooth and an oscillating 
part: 



E(N) = E(N)+E(N). 



(1) 



The smooth part describes the overall change in energy 
as the cluster size increases and is given by a liquid drop 
expansionoa 

E(N) = a x N + a 2 N 2/3 + a 3 N 1/3 + ■ ■ ■ . (2) 
The oscillating part is responsible for the shell structure. 



The idea of shell correction methods is to give a pre- 
scription for determining E(N) from a one-particle cal- 
culation. These methods were pioneered by Strutinsky, 
who showed how the oscillating part of the total energy 
resulting from a Hartree-Fock calculation for atomic nu- 
clei can be determined from the sum of the sioftLe^particle 
energies ^ e M of a suitably defined potentialO'tS A sim- 
ilar result holds for ELN) extracted from local-density 
functional calculatipBsJLj The latter are more common 
for metal clusters Jj'EJEj For clarity and to fix the nota- 
tion we give a short outline of the relevant argument. 

To find the ground state energy of a system of N elec- 
trons using densiUi-iunctional theory we use the Kohn- 
Sham formalism E3'Ej Starting from some electron density 
no(r) we have to solve the Kohn-Sham equations with the 
potential 



UmTS : i;,,(Fi + / d 3 r'^Q- + V xc [n ] 



(3) 



Having found the N lowest eigenstates ipp,{f) with ener- 
gies e M , an estimate of the total energy of the system is 
given by the variational expression 

E[n ] = E km [no] + E Co ui[no] + E xc [n ], (4) 
where the kinetic energy is given by 



JV 



(5) 



d 3 r V xc [n ] n(r ) - / d 3 r V ext (f) n(r) 



and the Coulomb energy is the sum of the Hartree energy, 
the interaction of the electron density with the external 
potential (e.g. the potential arising for the ion-cores), and 
the electrostatic self-energy of the ionic cores 



Ecoui [no] 



1 



d 3 rd 3 r' "ff"^ 



+ / (TrVestWnW + Ei. 



(6) 



The 'new' electron density n(r) in the above expressions 
is given by IVv(j)| 2 - 

Let us assume that n (r) was chosen close to self- 
consistency. Then n(f) will not differ too much from 
n Q (f) 



n(r) = n (r) + dn(f), 



(7) 



and we can expand the expression for the total energy 
(Q) in powers of Sn. Using 

1 I d 3 rd 3 r> r^r = ~l ld 3 rd 3 r> 



d 3 rd 3 -' 



\r — r | 



and 



2 



E xc [n] = E xc [no] + J 



,3 SE xc [n ] 



6n(f) + 2 (5n) 



we can, to first order in 5n, write the total energy as a 
functional of only the initial electron density no(r): 

^ ^ 9 / r — r> 



d 3 r V xc [no](r) n (f) + E xc [n ] + E T . (8) 



A good choice for n is the electron density utf re- 
sulting from an extended Thomas-Fermi (ETF) calcula- 
tion. Since n,ETF{N;r) varies smoothly as the number 
of the valence electrons N in the cluster is changed, all 
terms in (||), except for the first one, contribute exclu- 
sively to the smooth part E(N) of the total energy. I.e. 
to first order in Sn, all electronic shell effects E(N) are 
contained in the sum of the one-particle energies X) e M- 
Hence the oscillating part E(N) of the total energy can 
be determined from the spectrum of the family V(N; r) 
of Kohn-Sham potentials which arise from the electron 
density tijstf {f) ■ More generally, the above reasoning 
holds for all families of electron densities no(N;r) that 
are close to self-consistency and smooth in N . 

It is common practice to immediately work with pa- 
rameterized potentials V(N;r). Usually they are cho- 
sen to fit cxaapinxints or the results of self-consistent 
calculations .ouO'Ej Imagining that these potentials 
arise from a hypothetical family of electron densities, 
the above arguments still apply. The prototype of 
such a phcnamenological shell model is the Woods-Saxon 
potentialEM 



V(N;r) = 



-Vo 



1 + cxp ((r - r s N 1 / 3 )/^ ' 



Variants are the Wine-bottle potential!! a»d the Woods- 
Saxon potential with asymmetric surfacelilJ 



V(N;r) = 



-V 



+ exp ((r - (r s + AR))/f(r)) ' 



(10) 



where f(r) is an analytical function modeling the poten- 
tial near the cluster surface. 

To get a feeling for the approximations involved, we 
compare the results of a self-consistent calculation for 
gallium clusters to the oscillating part of the total en- 
ergy found using a family of model potentials (Fig. [I]). 
For the self-consistent calculations we used the homoge- 
neous, spherical jellium model. LI The potentials for the 
one-particle calculation were obtained by fitting a func- 
tion of the type JlC| ) 



V(N;r) 



-V 



l + exp((r- R(N))/a(r)) ' 



(11) 




JLLI 



FIG. 1. Comparison of the oscillating part E(N) of the 
total energy obtained from a self-consistent calculation using 
the homogeneous, spherical jellium model (dotted line) to the 
E(N) extracted from the sum of the one-particle energies of 
a family of model potentials (full line). The model potentials 
are given in equation (|Io|), with parameters r s — 2.19 ao, 
Vo = 1.04 Ry, AR — 0.73ao, ao = 1.03 ao, ai = 1.13 ao, and 
a2 = 0.21 ao. 
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FIG. 2. Oscillating part of the total energy and of con- 
tributions to it (cf. eqn. (Q)), obtained from a self-consistent 
calculation using the homogeneous, spherical jellium model. 
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where R(N) = r s N 1 / 3 + A(r) and a(r) = a + 
ai tanh(a2(r — R(N))), to the self-consistent Kohn-Sham 
potentials for jellium clusters having 1500, 3000, 4500, 
and 6000 valence electrons. 

To emphasize the importance of the smoothness of the 
model potentials V(N; r) in TV, and to demonstrate from 
what subtle cancellations the electronic shell structure 
arises in self-consistent calculations, we show in Fig. || 
the oscillating part of the total, the kinetic, the Coulomb, 
and the exchange-correlation energy determined from a 
jellium calculation for gallium clusters. Even though in 
this calculation we are dealing with the self-consistent 
potentials and electron densities (i.e. 8n = 0), Eki n {N), 
which contains the sum of the Kohn-Sham energies, is 
by far not the only term contributing to E(N). More 
surprisingly the electronic shell structure as revealed by 
E cannot be found in any single contributions to the to- 
tal energy. For the self-consistent calculation it rather is 
resulting from the subtle interplay of the different oscil- 
lating terms. 



III. PERIODIC ORBIT EXPANSION (POE) 

As we have seen in the previous section, the oscillat- 
ing part E of the total energy can be extracted from the 
sum £ i °f the N lowest eigenenergies for a suitably 
chosen family V(N; r) of model-potentials. The deter- 
mination of the electronic shell and supershell structure 
is thus reduced to an eigenvalue problem for these po- 
tentials. Furthermore, the radius of the clusters we are 
interested in is considerably larger than the de Broglie 
wave-length of the electrons at the Fermi level. The 
semiclassical approximation seems therefore well suited 
for solving the single-electron problem in question. In 
fact, for the spherical cavity, a simple rescaling of the 
Schrodinger equation shows that the limit R — > oo is 
identical to the semiclassical limit % — > 0. 

The salient feature of the semiclassical approach to de- 
termining the electronic shell and supershell structure 
is, that it provides a natural splitting of the density of 
states (and consequently the total energy) into a smooth 
and an oscillating part. The smooth part corresponds to 
Thomas-Fermi theory, while the quantum corrections are 
given by a sum over the nontrivial periodic orbits. For 
an understanding of the oscillating part E of the total 
energy, we need only consider the latter. 

In the present section we derive the periodic-orbit- 
expansion (POE) for spherical potential wells. Starting 
from the oscillating part p(k) dk of the density of states 
for a given potential V(r), we proceed to a POE for the 
oscillating part E(N) of the total energy for a family 
V(N; r) of potentials. We illustrate the results by giving 
explicit expressions for infinite potential wells. These will 
be the point of reference for the leptodermous expansion 
discussed in Sec. [v|. We furthermore use the spherical 
cavity to assess the validity of the various approximations 



made, by comparing the semiclassical results to the re- 
sults obtained from numerically solving the Schrodinger 
equation. 

A. POE for the density of states 

We first sketch the derivation of the periodic orbit ex- 
pansion of the density of states for a spherical poten- 
tial well al©ngi th* lines of the path-integral approach of 
GutzwillerE'EirEil The starting point of the derivation is 
the relation between the density of states p(E) and the 
Green's function G(E) 

p(E)dE = — 3fTr G{E +ie) dE. (12) 

Expanding the trace in real space leads to 

p(E)dE = —9 / d 3 r lim G(f,f ;E + ie)dE. (13) 

Using the transformation 

7 f°° 

G(r,r Q ;E) = --J^ dte iEt ' n G(f,f ;t) (14) 

the energy-dependent Green's function is expressed in 
terms of the time-dependent Green's function, which in 
turn can be written as a path-integral 

G(f, Fo; t) = J V[t{t)] e lS ^' h . (15) 

Here the integration is over all paths r(r) in configura- 
tion space connecting the point fo with f taking a time 
t. S[f(r)] is the classical action along such a path. In the 
semiclassical limit Ti — > only those paths contribute for 
which the phase is stationary (SS — 0: classical paths). 
If the second variation of the classical action is finite 
(5 2 S 7^ 0), the path-integral becomes a Gaussian inte- 
gral which can be evaluated analytically. For S 2 S = 
the path-integral picks up a phase factor. Thus in the 
semiclassical limit the pathjptegral (|l5|) decomposes into 
a sum over classical pathsoO 

Mm G{r,r ;t) = (16) 

H— >0 

The transformation ( |l4| ) of the time-dependent into the 
energy-dependent Green's function can also be evaluated 
by the method of stationary phase, leading to an expres- 
sion of the form 

G(f,f ;E)= J2 ^ exp fiS[f] + <ty . (17) 

Again there are additional phases associated with the 
vanishing of the second variation. 
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Taking the trace ([l3]) involves integrating over ro an d 
taking the limit and r — > ro . The integration over ro is 
again done in the semiclassical limit. Here the stationary- 
phase condition requires that the final moment equals the 
initial moment. The limit r — > f, D finally closes the or- 
bits. Thus, since the orbits return to r*o with the same 
momentum, they are closed in phase space, i.e. they are 
periodic. There are two distinct classes of such orbits. 
The first consists of only the direct path, the length of 
which vanishes as r — > rb- It consequently is local and 
gives rise to the Thomas-Fermi density of states p. The 
second class consists of periodic orbits of finite length. 
These non-local paths give rise to a quantum-correction 
to p. Hence in the semiclassical approximation the den- 
sity of states is given by the local Thomas-Fermi term 
with non-local corrections described by a sum over peri- 
odic orbits: 

p{E) dE = (p(E) + p{E)^j dE. (18) 

In a spherical potential well, i.e. a potential with at 
most two radial turning points, all periodic orbits can be 
easily enumerated: A periodic orbit is characterized by 
the number A of times it winds around the origin and the 
number v of times it traverses the outer turning point. 
By symmetry all orbits (A-i/) that only differ in orien- 
tation are equivalent. Fig. |3| shows some of the periodic 
orbits for a spherical cavity. The periodic orbits for a 
general spherical potential well are mere rounded, but 
are still described by the pairs (A, z/).ErE3 

The periodic orbit expansion for the oscillating part 
of the density of states is thus given by a sum over the 
families of equivalent orbits (A, v) 




(1,2) (1,3) (1,4) (1,5) (1,6) (1,7) 




(2,4) (2,5) (2,6) (2,7) 




(3,6) (3,7) 



FIG. 3. Some periodic orbits for a spherical cavity. They 
are characterized by the pair (A, v), where A denotes the 
number of times the orbit revolves around the origin before 
it closes on itself, and v is the number of vertices it has. 
Note that (nX, nv) is the orbit obtained by traversing (A, v) 
n times. 

the phase is given by 

<P(\,v) = (%> + A - ~) 7T, (22) 
and the amplitude takes the form 

A(\,u) = VkRl /2 ot(x,u) (23) 
with the dimensionless geometry-factors 




p(E) dE 



S, 



ip {x ,u) dE, (19) 



where Sf\,u) is the classical action for an orbit (A, v) and 
f{X,v) is the Maslov phase. The amplitude with which 
the orbit (A, v) contributes is given by 



(X,v) 



717/ 



ds r /h 



d 2 s r h 



dL 2 



-1/2 



(20) 



with L-denoting the angular momentum and s r the radial 
action Jj 

For a spherical cavity of radius Rq the terms that en- 
ter the periodic orbit expansion take a simple form: The 
classical action of an orbit equals its length times the 
wavevector k 



Si\ jl/ )(k)/h = 2i/kRo sin 



ttA 



(21) 
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FIG. 4. Amplitudes Q(a,i/) with which the periodic orbits 
(A, v) contribute to the oscillating part p of the density of 
states. Shown are the amplitudes (crosses) for no (fo = 0) 
and for an intermediate (fej = 0.1) smoothing. To guide the 
eye, the amplitudes for a given number of turns are connected 
by lines. 
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The cx(\ t v) determine the relative importance of the 
periodic orbits in the POE ([jl]). Their values for the 
first few periodic orbits are shown in Fig. ^ (fej = 0). 
We note that the amplitudes for the linear orbits (A, 2A) 
vanish. This can be understood by a simple dimensional 
argument .Q Since the sum in the POE is over all periodic 
orbits, the number of different but equivalent orbits (A, v) 
will be reflected in the amplitude ot(\ tV \. To parameterize 
all the different orientations of the linear orbits it is suffi- 
cient to give the coordinates of one of their outer turning 
points. Since for a spherical potential well the outer turn- 
ing point lies on the surface of a sphere, the manifold of 
the linear orbits has dimension 2. All the higher orbits 
are not linear but lie in a plane, so we need an addi- 
tional parameter to fix the orientation of this plane. The 
manifold of the planar orbits are therefore 3-dimensional 
('there are many more planar than linear orbits'). Thus 
the linear orbits do not contribute to the leading order of 
the periodic orbit expansion. The largest amplitudes are 
found for the triangular (1, 3) and the square (1, 4) orbit. 
The contribution from other orbits is, however, still large, 
i.e. one has to include many periodic orbits in a partial 
summation of (|f^) before one obtains a result close to p. 
This slow convergence is to be expected since the den- 
sity of states for a finite system is given by a sum of 
(^-functions which cannot easily be reproduced by a sum 
of analytical functions. To improve the convergence of 
the expansion one can replace the (5-peaks in the DOS by 
Lorentzians of width 7. This corresponds to introducing 
a complex wavevector k = k r + ik t in the periodic orbit 
expansion. As can be seen from Fig. [|, a finite value of k% 
serves to reduce the contribution of higher orbits consid- 
erably. However, since the shell and supershell structure 
in metal clusters is not directly linked to p but rather to 
the variations E in the total energy, we proceed to derive 
a periodic orbit expansion for E. As we will see, such an 
expansion converges much more rapidly than that for the 
density of states. We therefore need not introduce any 
smoothing. 



B. POE for the total energy 

To find a periodic orbit expansion for the oscillating 
part E of the total energy using the POE for p, we start 
from the integral 



<-E F (N) 

E(N) = I Ep{N; E) dE, 



(25) 



where the Fermi energy Ep(N) is fixed by the number of 
electrons N in the cluster 



pEp(N) 

N — p(N\E)dE. 
Jo 



(26) 



Similar equations hold in Thomas-Fermi theory. Sub- 
tracting the corresponding Thomas-Fermi expression 
from d2q) we find 



= / p(N;E)dE+ / p(N;E)dE. (27) 

Je f (n) Jo 

Since the smooth part of the density of states does not 
vary much over the small interval Ep . . . Ep, we can ap- 
proximate the first integral in the above expression by 
Ep(N) p(N; Ep). We then can use the above equation to 
solve for Ep(N). In a similar fashion we can approximate 
the difference of ( ^5| ) and its Thomas-Fermi counterpart 
by 



E F (N) 



E F (N) 



Ep(N; E) dE. 



E(N) re Ep(N)Ep(N)p{N-Ef) 



Using the approximate expression for Ep from eqn. (|2 
we find 



E(N) 



E F (N) 



(E - E F {N))p{N-E)dE. 



(28) 



Since the integrand vanishes at the upper limit of inte- 
gration it is now possible to approximate Ep(N) by its 
Thomas- Fermi counterpart Ep(N). Integrating by parts 
we finally arrive at 

pE F (N) j-E 

E(N) pa — / dE dE' p{N-E'). (29) 

Jo Jo 

i.e. to find an approximation to the oscillating part of the 
total energy we have to integrate twice over the oscillat- 
ing part of the density of states. 
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FIG. 5. Amplitudes f3(\, v ) with which the periodic or- 
bits (A, v) contribute to the oscillating part E of the total 
energy. Comparison with Fig. ^ shows that the dominance 
of the short, planar orbits, which for the oscillating part p 
of the density of states has to be enforced by introducing an 
artificial smoothing ki, occurs naturally for E 
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5 10 15 20 25 30 

N 1/3 

FIG. 6. Comparison of the periodic orbit expansion for the oscillating part E of the total energy to the quantum mechanical 
result for spherical cavities. The plots show E obtained from a truncated periodic orbit expansion, including more and more 
orbits. The orbits were included in the order of decreasing amplitude: (1,3), (1,4), (1,5), (1,6), (1,7), (1,8), (1,9), (2,6), 
(2, 5), (2, 7), (1, 10), (2, 8), (1, 11). The Eqm from the quantum mechanical calculation is shown by the dotted line. 



Using ( p9| ) and ( |l9| ) we find for spherical cavities of 
radius Ro(N) the expansion 



E(N) 
\/k F R k% 



(30) 



4a 



S? 



(A, i/) (A,i/) 



) k F R - tp(x,v) 



which is similar to (|19[), the main difference being the 
change in the amplitudes: Due to the twofold integration 
the amplitudes are divided by the square of the dimen- 
sionless classical action Sr\ jU \ = Sn v -\/(KkRo). The new 
geometry factors are thus given by 



4a 



o2 



(31) 



They are plotted in Fig. |[ A comparison with the cx(\, v ) 



(Fig. ||) shows how the contributions of the long orbits 
(with large classical action, see eqn. (pi])) to the POE for 
E(N) are reduced. This improvement of convergence can 
be understood intuitively since E(N) is continuous, while 
the density of states is highly singular, being a forest of 
((-functions. 

To check the approximations made in the derivation 
of eqn. (p9| ) we compare the results of a truncated peri- 
odic orbit expansion for E with the oscillating part Eqm 
of the total energy derived from a quantum mechanical 
calculation. Such a comparison for spherical cavities of 
radius Rq = N 1 / 3 is shown in Fig. ^. It turns out that 
the truncated POEs reproduce Eqm very well, even if 
only a few periodic orbits are included. In particular the 
first two supershells can be described using only the tri- 
angular and square orbit. However, for even larger sizes 



7 



N it seems that higher orbits are needed to describe the 
structure in Eqm- We note that due to its nature of 
being a semiclassical result, the periodic orbit expansion 
for E will not converge to Eqm but to its semiclassical 
approximation. 



IV. LATTICE CONTRACTION 



an isolated factor k F . Since the Fermi energy obvi- 
ously depends on the cluster size, we expect the ampli- 
tude of E to be affected by a lattice contraction. For 
R(N) = r s A 1 / 3 + AR we find from (||) 



k F (N)r s =61+ 62 



bl —) +0(N- 2 ^). 
r s 



(35) 



As a first application of the periodic orbit expansion 
for the oscillating part of the total energy we look at 
the effects of a lattice contraction on the electronic shells 
and supershells. Extended X-ray absorption fine struc- 
ture (EXAFS) measurements on small clusters adsorbed 
on some substrate indicate that the next-neighbor, di g - 
tance might decrease with decreasing cluster sizeO'EHl 
Such a contraction of the cluster is also suggested by 
a continuum description of large clusters: Since for small 
clusters the surface-to- volume-ratio increases, the surface 
tension, which has the tendency to compress the cluster, 
will become more and more important. Neglecting this 
finite-size effect, a spherical cluster of iV atoms will have a 
radius Rq = rsA 1 / 3 , where r s is the Wigner-Seitz radius 
of the bulk material. Acting against the compressibility 
k, the surface tension a will reduce the cluster radius: 



AR _ 2na 1 
Ro 3 Rq 



(32) 



Thus we expect the radius of a cluster of A atoms to be 
given by 



R(N) = R a + AR = r s A 1/3 + AR 



(33) 



with a size-independent contraction AR. A rough esti- 
mate using the compressibilities and surface tensions of 
the liquid alkali metals just above the melting point gives 
Ai?/r s » 0.1. 

To see how the electronic shells and supershells are af- 
fected by such a lattice contraction we first turn to the 
special case of the spherical cavity. Inspecting eqn. j30| ) 
we see that the periods of the oscillations in E are de- 
termined by the product k F R, while all other quantities 
in the argument of the cosine are independent of clus- 
ter size. k F R should, however, depend only weakly on 
the cluster radius, since we expect a decrease in cluster 
radius to be compensated by an increase in the Fermi 
energy. For the spherical cavity we actually have 



k F R = 61A 1 / 3 + 6 2 + C(A -1 / 3 ) 



(34) 



with bi = (9-71-/4) ^ and 6 2 = 3ir/8. I.e. for spherical 
cavities k F (N)R(N) is independent of R(N) and there- 
fore the position of the electronic shells and supershells 
are independent of the actual cluster radius R(N). 

Looking again at eqn. (^) we notice that in the over- 
all prefactor of the periodic orbit expansion we have 



Thus for small cluster sizes the amplitude of the shell os- 
cillations will increase with increasing lattice contraction. 

To check in how far the above findings also hold for 
more realistic cluster potentials, we have calculated the 
oscillating part of the total energy from the eigenval- 
ues of a family of Woods-Saxon potentials. The param- 
eters were chosen to resemble the potentials from jel- 
lium calculations for sodium. In Fig. |?] we compare Eqm 
for contracted/expanded clusters with clusters of radius 
R(N) = r s A 1 / 3 . We find that even for unphysically large 
contractions/expansions (AR/r s — ±0.5) the location of 
the electronic shells and supershells is hardly changed, 
while for small numbers of atoms A the amplitude of 
E(N) increases/decreases for the contracted/expanded 
clusters. 




FIG. 7. Oscillating part E of the total energy for 
Woods-Saxon potentials V(r) = —Vo/(l+exp((r — R(N))/a) 
with lattice contraction R(N) = rsN 1 ^ 3 + AR. The parame- 
ters r a = 3.93 ao, Vo = 0.46 Ry, and a = 0.94 00 were chosen 
to resemble the potential for sodium-jellium clusters. 
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V. LEPTODERMOUS EXPANSION 

So far we only have explicit expressions of the periodic 
orbit expansion for cavity potentials. We now want to 
extend the POE to more realistic potentials V(r), like 
the Woods-Saxon potential (eqn. (||)). These potentials 
differ from the cavity potential by having a surface of fi- 
nite width a. Since the slope of the cavity potential is 
infinite at the surface, SV(r) = V(r) — V ca vity is never 
a small quantity. But, rewriting integrals over V(r) in 
a suitable way, we can use the surface width a as an 
expansion parameter. 

From (^) we see that we need to find expressions for 
(i) the action integral «5(a,j/) for orbits (A, v), (ii) the 
Fermi wave-vector Uf in Thomas-Fermi approximation, 
and (iii) the phases <P(a,i/) • For this we proceed as follows. 
We first introduce the idea of the leptodermous expansion 
for integrals over the potential V(r) for classical action. 
We find an expansion 



S 



qcavity , r, f 
(A, i/) + iVli 



a 

Rq 



(36) 



Then we estimate the change in kpr s due to the finite 
surface width. To first order in a we find 



kvT. 



Finally we estimate the phase y>(\ tV )- Rearranging terms 
in powers of TV 1 / 3 , the argument of the cosine in eqn. 
( p0| ) then reads 



9tt 



1/3 



§c*v«V N l/3 

4 / (A,w) 



(38) 



+S, 



cavity 
(A,*) 



Cl + c 2 - 



2v 



n \ 1/3 

9tt\ ' ~ a 



i.e. the first order terms in the leptodermous expansion 
give rise to a phase shift in the periodic orbit expansion, 
while the frequencies S?^ v \ y are unchanged. 



A. Classical action 

To introduce the basic idea of the leptodermous ex- 
pansion we first consider potentials that differ from the 
cavity potential only in a small region around the cluster 
surface, say for r > Rq — a. It is then straightforward 
to split the radial integrals into two parts, one integral 
over the interior r — . . . Rq — a and one over the surface 
region r = Rq — a . . . r out . Thus the radial action can be 
rewritten as 



s r /h 



^E-V(r) - L 2 /r 2 dr 



dr + yC dr, 

J Ro—a 



(39) 



The first integral is the action integral for a spherical cav- 
ity of radius Rq — a, which we know already. The second 
integral can, in general, not be solved analytically. But 
for small a we can get a good approximation by neglect- 
ing the variation of the angular momentum term over the 
small interval [Rq — a, r out ], e.g. setting L 2 /r 2 to L 2 /Rq. 

Realistic cluster potentials are not that simple. In a 
Woods-Saxon potential there is no obvious point, that 
separates bulk from surface. We can still make the 
same ansatz by choosing some small a, but we have to 
make sure that our result does not depend on our spe- 
cific choice. To do so, we add and subtract the integral 
J* " ^E + V Q ~L 2 /R 2 dr to ©. Using V(r) « -V 
for r < Rq — a we find 



s r /h 



Ro-a 







Ro-a 



yJE + V - L 2 /r 2 dr 



E-V(r) - L 2 /R 2 dr 



E + Vq- L 2 /R 2 dr 



(40) 



Using again L 2 /r 2 w L 2 /Rq for r e [R — a,r out ], we 
can extend the upper limit of integration of the first and 
third integral from Rq — a to r out . The first integral is 
then the radial action for the cavity potential, the two 
other terms give the correction due to the soft surface 
of the potential V(r). Introducing dimensionless quan- 
tities a = a/i?o, s — s/fci?o, and L = L/hkRo, with 
k = \JE + Vq, and expanding in powers of the reduced 
surface- width a, we get 

s r (P, L, a) = s c r amty (L) + t(P, L)a + 0{a 2 ), (41) 



where P = ^{E + Vq)/Vq, and s^. avity (L) = \J\-L 2 - 
Larccos(L) is the reduced radial action for a spherical 
cavity. In appendix ^ it is shown how to calculate / s for 
a Woods-Saxon potential. 

Given the expansion for the radial action, we now pro- 
ceed to calculate the action for a periodic orbit (A, v) 



2v s r + 27rA L 



(42) 



The angular momentum Lr\ tI/ ) associated with the orbit 
can be determined from the periodicity condition: In or- 
der to close after A turns and having traversed the outer 
turning point v times, the angle $ swept during one ra- 
dial oscillation must be tcX/u. With (Bfi) this leads to 



— = $ = 
v 



^ = arccos(f) -^±a + 0(a 2 ). (43) 



dL 



Taking the derivative with respect to a at a = 0, we can 
solve for the first order correction in the reduced angular 
momentum: 
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L 



r cavity 



, f cavity^ 



di s {p,L c ( a ^y) 



i)L 



(44) 



leptodermous expansion is exact, i.e. there are no higher 
order terms in ({4l|). From (|47| ) the surface parameter 
thus is 



We can use this to expand s c r V (L) in (41) around 



L C {\ % v) ~ cos{n\/v). Inserting into ( |42| ) we see that the 
first order correction in §^ amty cancels that coming from 
L. Thus to first order in a the reduced action for a peri- 
odic orbit (A, v) is given by 

S(P,L {x>v) ,a) = S™«» + 2 V i s (P,L™^)a + 0(a 2 ). 

(45) 

This result is independent of the specific form of the po- 
tential, as long as an expansion (|4l]) of the radial action 
exists. 



c F 



3a 



I s (P(k),0) k 2 dk 



1/3 



-It)' 



(52) 
(53) 



The analytic expression of In for Wo ods-S axon poten- 
tials is calculated in appendix [X| (eqn. ( A10 )). 



It is interesting to note that we can obtain the same re- 
sult from expanding the Thomas-Fermi integral in powers 
of a 



B. Fermi level 

We now turn to the problem of determining the Fermi 
wavevector kp in the extended Thomas- Fermi approx- 
imation for a cluster with N electrons. In an infinite 
system we have k F r s = (9-7T/4) 1 / 3 . For a finite system 
there will be corrections arising form the surface 

k F r s = (9tt/4) 1/3 + c f N~ 1/3 + ■ ■ ■ (46) 

The surface term can be calculated from the quan- 
tum mechanical scattering phase ip(k) at the surface 
potentialEZru 

CF = ~W Jo" (J ( 4? ) 

Assuming that the potential is slowly varying at the outer 
turning point, we can determine the scattering phase 
from the classical action. 

For a slowly varying potential the WKB wave function 
in the region, r Rq — a where the potential is practically 
constant isEj 

uwkb(t) oc cos ^ k(r) dr - , (48) 

while the quantum mechanical wave function is 

uqm{t) oc cos(fc(iio - r) - <p{k)). (49) 

In the semiclassical limit, i.e. for large Ro, both expres- 
sions should be equal. Choosing r = 0, we find 

tp{h) = -\J ^ fc M dr - k R o) + \ ( 5 °) 
= j~I s (P,0)ka + O{a 2 ), (51) 

where in the last equation we have used the linearization 
([ll]) of the radial action. We note that for L = the 



N= (E F -V{R)f /2 r 2 dr 

Jo 

= { -^f^ + (k F R ) 3 I N a (54) 



and solving for kpr s . In fact, in the approximation con- 
sidered here, both approaches are equivalent. 

It is important to realize that the above reasoning in- 
volves two, possibly conflicting approximations. In the 
ansatz (^8|) for the semiclassical wave function we have 
assumed that the potential is slowly varying, i.e. that a 
is large enough, while for the leptodermous expansion of 
the radial action of the orbits with L ^ we require that 
a is small. 

To see how well the expression ( f52| ) for cp works, we 
compare it to the surface term obtained by fitting the 
Fermi wavevector kp(N) calculated quantum mechani- 
cally for Woods-Saxon potentials holding 50 . . . 8000 elec- 
trons. As expected our approximation approaches the 
quantum mechanical result in the limit of large surface 
width a (slowly- varying-potential regime). But it also 
works quite well for relatively small a. For very small 
a the approximation of course breaks down, since our 
ansatz does not describe the crossover from the slowly- 
varying-potential regime to the potential step at a = 0, 
for which the phase in uwkb{t) is arctan(«;/fc) instead of 
7r/4. Using this phase in the above derivation, we recover 
the correct surface parameter for the finite potential well. 

Since for larger a the approximation to c F runs roughly 
parallel to the true value of the surface term, on might 
improve the accuracy of the leptodermous expansion by 
shifting (|5^) by a constant (a-independent) amount c\. 
Since the Woods-Saxon potential in one dimension is ex- 
actly solvableE3e3, this is straightforward. 
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FIG. 8. Surface term cf in the asymptotic expansion 
of kpr s (cf. eqn. (^)) as a function of the surface width a 
for a Woods-Saxon potential. The potential parameters are 
Vo = 0.45 Ry, r s = 4 ao)- The crosses were obtained by fitting 
the Fermi wavevector from full quantum mechanical calcula- 
tions for clusters with 50 . . . 8000 electrons. The dashed line 
gives the estimate of cf as given by eqn. (p2|). 



the periodic orbit expansion ( pCf ) of the oscillating part E 
of the total energy. The frequencies associated with the 
periodic orbits turn out to be unchanged, to first order 
there is only a phase shift: 



E(N) cx £ /J (A> „, cos f ffl j + A* (Afl0 

(57) 



1/3 



with 



A$ 



9?r 

T 

W + \ 



1/3 



2^(P,L)-S ( ^/ W (P) 
1/4]tt. 



a 
r. 



(58) 



C. Maslov phase 

For separable potentials the Maslov phase is deter- 
mined by the phases that the semiclassical wave func- 
tion picks up at the classical turning points. For a given 
periodic orbit (A, v) there are 2A turning points in the 
^-motion, each contributing a phase tt/2. The same or- 
bit also has 2v radial turning points. The v inner turn- 
ing points see the smooth centrifugal potential, hence 
also contribute tt/2. The phase <j) ut at the outer turning 
point depends on the shape of the potential V(r) at the 
surface. The Maslov phase for the orbit is then given by 

<P(\, V ) = [(1/2 + <j>out/K)v + A — l/4]vr. (55) 

For a step potential, we can find (pout by matching 
the semiclassical radial wave function to the boundary 
condition at the turning point. For an infinite poten- 
tial well (pout = Tj while for a well of depth Vq (pout = 
2 arctan(«;/fc), where k — \JV$ — k 2 . 

For a slowly varying potential, on the other hand, 
the standard result obtained by linearizing the potential 
around the classical turning point is (pout = tt/2. 

The Woods-Saxon potential for a typical cluster has 
a surface width of a = 0.5 . . . 1.5 ao. As we have seen 
above, for such values of a the potential is already in the 
slowly- varying-regime. We therefore use 

f(x.u) ={v + X- 1/4]tt. (56) 



D. Leptodermous POE 

We can now collect all the contributions to calculate 
the effect of a softening of the potential at the surface on 



We stress again that we have made two, possibly con- 
flicting approximations. On the one hand, the leptoder- 
mous expansion of the radial action relies on the fact 
that the surface width a is small, while a slowly- varying- 
potential assumption enters in the calculation of the scat- 
tering phase. To see how the above expression works 
in practical calculations, we compare the periodic orbit 
expansion (|57j) with the result of quantum mechanical 
calculations for Woods-Saxon potentials with parameters 
typical for alkali metal cluster, see Fig. ^|. The agreement 
is surprisingly good. The shift of the supernodes with 
increasing surface width is well described, and also the 
shell oscillations are quite well reproduced. We could get 
even better agreement by numerically fitting the action 
integrals (see Fig. [ll]) and the surface coefficient cf with 
linear functions in a. In that sense the concept of the sur- 
face introducing just a phase-shift in the periodic orbit 
expansion seems to be applicable even beyond the range 
where the analytical expressions from the leptodermous 
expansion are good approximations. 

Taking only triangular and square orbits into account, 
the shift in the shell (supershell) oscillations is given by 
(A$(x,3) =t ^ < J > (i.4))/2- As we can see from (j5q), the con- 
tributions coming from cf almost cancel for the super- 
shells, since the classical actions for the triangular and 
the square orbit are so similar. Because of this can- 
cellation of errors the shift of the supershells with the 
surface width is very well described in the leptodermous 
expansion.ES The shell oscillations are more sensitive to 
approximations. But the most important feature, namely 
that shells are hardly affected by changes in the surface 
width, is also well reproduced. 
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5 10 15 20 

N 1/3 

FIG. 9. Oscillating part E of the total energy for 
Woods-Saxon potentials V(r) = -V /(l + exp((r- R(N))/a) 
with different surface width a. R{N) = N 1 ^ 3 r a with r 3 = 4ao 
and Vo = 0.45. The dotted lines give the results of quantum 
mechanical calculations. The full lines are obtained from the 
leptodermous expansion ([37|), including only triangular and 
square orbits. 



VI. CONCLUSIONS 

Using periodic orbit theory we see that the electronic 
supershell structure is a sensitive probe for the surface of 
clusters. It turns out that possible lattice contractions, 
which at first sight seem like important surface effects, 
hardly influence the electronic shell structure. There is 
a pronounced effect of the width of the surface region 
on the position of the supershclls. The leptodermous ex- 
pansion around the limiting case of a spherical cavity 
provides a natural framework for understanding the shift 
of the supershells with increasing surface width. A par- 
ticularly nice feature of the leptodermous expansion, as 
we have presented it here, is the fact that the expan- 
sion of the radial action is the only input we need. All 
other quantities entering the periodic orbit expansion can 
be easily derived from the radial action. It is therefore 
straightforward to apply the formalism to other types of 
potentials. 

The shift in the electronic supershells that is described 
by the leptodermous expansion has hcpn seen in numer- 
ical studies of E for soft potentialaEl and it has been 
used to understand the results of self-consistent jellium 
calculations .Ej The observations of a shift proportional 
to the surface width can be regarded as a signature of 
the leptodermous regime. Eventually the leptodermous 
approximation will break down, since for extremely soft 
potentials the planar orbits with A = 1 cease to exist. 



For such potentials star orbits become the leading terms, 
which causes a changp-iia, the frequency of the shell and 
supershell oscillationscIle3 as opposed to a change in the 
phase only. 

It is interesting to compare the leptodermous expan- 
sion of the semiclassical sum over periodic orbits with 
the quantum mechanical perturbation theory. In quan- 
tum mechanics we would expect perturbation theory to 
break down when the shifts of the energy levels are of 
the order of their spacing. For the potentials we have 
considered here the change in the energy levels is quite 
large, especially for the levels with high angular momen- 
tum, which are most sensitive to the potential at the sur- 
face. Typical shifts of the energy levels with the surface 
width a for a set of Woods-Saxon potentials are shown in 
Fig. [HJ Nevertheless, the electronic shells are not that 
strongly affected, because the levels with large angular 
momentum, which are mostly responsible for the oscilla- 
tions in the total energy, are shifted by large, but similar 
amounts. It seems that since the semiclassical periodic 
orbit expansion does not deal with individual energy lev- 
els but only with the collective changes in the spectrum, 
it works so well, even for large surface widths a. 




FIG. 10. Shift in the energy levels e n ,i for Woods-Saxon 
potentials with different surface-width a. The parame- 
ters for the potentials are Vo = 0.45 Ry, Ro = 50 oo and 
a = 0.5 . . . 1.5 ao. For given angular momentum I the levels 
are plotted one above the other. Each line shows e n j as a 
function of a, with a increasing from left to right. 

The leptodermous expansion should also be of use in 
understanding experiments probing the transport prop=. 
erties of high- mobility semiconductor microstructures.EJ 
The oscillations in the conductance of quantum dots are 
quite well described by simple cavity potentials. At first 
sight this seems surprising since the confining potential of 
a quantum dot is expected to be rather smooth. However, 
if these potentials are still in the leptodermous regime, 
it is clear that calculations using simple cavity potentials 
(or billiards) already describe the qualitatively correct 
physics. 
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APPENDIX A: LEPTODERMOUS EXPANSION 
FOR A WOODS-SAXON POTENTIAL 

In this appendix we show how to evaluate the inte- 
grals in the leptodermous expansion for a Woods-Saxon 
potential (Q). 



arctany: = arcsin ^.Py 1-L 2 j + 0{a 2 ). (A5) 

The first term is a bit more difficult, since for a — > the 
artanh^/7 diverges as l/2a. We find 



2aartanh^/7 = 1 + 2d In \^2P^ l-L 2 j + 0{d 2 ). (A6) 

Collecting our results we see that we can expand the 
radial action for a Woods-Saxon potential with surface 
parameter a around the radial action of the correspond- 
ing spherical cavity. The term of first order in a is given 
by I S (P, L) a with 



1. Radial action 

For a potential with small surface width a the radial 
action can be approximately written as (see discussion of 
eqn. ©) 



S r /h : 



y/E + V a - L 2 /r 2 dr 

Tout 



(Al) 



E -V(r) - L 2 /R 2 dr 



Ro 



E + Vq - L 2 /R 2 dr. 



The first integral is the radial action for the spherical 
cavity. The third integral is trivial: the integrand is a 
constant. Introducing k = y/E + Vb and L = L/hkRo to 
simplify the notation, it is given by hkRo Vl -L 2 . The 
second integral is more difficult. Rewriting the Woods- 
Saxon potential (P) as 



V{r) 



-Vb 



V 



1 + tanh 



Ro 



2 V \ 2a 

and substituting y = (r — Ro)/a, we are led to 



V2> 



artanh(c) 



\Jc- tanh(y) dy 



(A2) 



(A3) 



-l/2a 

with c = 2P 2 (1 — L 2 ) — I, P = ^/(E + V )/V , and 
a = a j ' Rq. This expression can be evaluated analytically 



2d I yl — L 2 artanh 



1 _j /2 _ l + tanh(-l/ a) 
2P 2 



1-L 2 



(A4) 



Try — (1 — L 2 ) arctan 



i 2 l + tanh(-l/2a) 



1-L 2 - 



1/P 2 -(1-L 2 ) 



Since the ansatz (eqn. (Al)) is already a first order ap- 
proximation, we need only the expansion of the above ex- 
pression for small a. For the second term this is straight- 
forward 



I s (P,L) = -[P L ln(2P L )- 



1 - P| arcsin(P L ) 



(A7) 



where we have introduced Pl = P^J 1 — L 2 . 



CO 




FIG. 11. Change of the radial action (in units of K) 
for Woods-Saxon potentials with different surface width a. 
The parameters for the potential are Vb = 0.45 Ry and 
Ro = 40 ao- The full lines give the radial action s r (E) for 
a triangular and a square orbit, k = (97T/4) 1 / 3 l/r 3 . The 
dotted line shows the change in the radial action for the tri- 
angular orbit as calculated using the leptodermous expansion 
(J s ka). The dashed line shows the same for the square orbit. 

To check how good this linearization of the radial ac- 
tion works in practice, we compare it to numerical results. 
Fig. [ll] shows such a comparison for a Woods-Saxon po- 
tential that roughly resembles a sodium cluster with 1000 
electrons. As can be seen, the leptodermous expansion 
works quite well. It is exact for L = (which is the 
quantity entering in the calculation of cp)- For the tri- 
angular orbit the approximation is very good even for 
rather large surface widths a. For the square orbit the 
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expansion still works well, although we start to see de- 
viations for larger a. This is due to the fact that in the 
leptodermous expansion we make an approximation in 
the angular momentum term, which becomes more crit- 
ical for orbits with large v. It is worth noting that the 
change in the radial action is surprisingly linear in a. 
That means, we could obtain results in the spirit of the 
leptodermous expansion (i.e. having only phase shifts in 
the periodic orbit expansion) , by fitting the results of nu- 
merical calculations of the action integrals with a linear 
function. Using such fits we could improve the accuracy 
of our results. 



2. Fermi level 



From equation (52) the surface term for kpr s in the 
leptodermous expansion is given by 



cf 



3 a 

^F JO 



/ s (P,0)fc 2 dk. 



(A8) 



Inserting (A7) we find 



cf = 



6Vo /2 a 



Pf 



P 2 ln(2P) - P\J\ - P 2 arcsin(P) 



dP, 



(A9) 

with P = k/^/Va- The first integral is straightfor- 
ward, the second is easily evaluated by substituting 
y — arcsin(P). We thus obtain eqn. (|53|) with 



I N (P) = 2 Mn(2P) 



p2 



.3/2 



arcsin(P) 
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